Bellwethers of change: population modelling of North Pacific humpback whales from 2002 through 2021 reveals shift from recovery to climate response

For the 40 years after the end of commercial whaling in 1976, humpback whale populations in the North Pacific Ocean exhibited a prolonged period of recovery. Using mark–recapture methods on the largest individual photo-identification dataset ever assembled for a cetacean, we estimated annual ocean-basin-wide abundance for the species from 2002 through 2021. Trends in annual estimates describe strong post-whaling era population recovery from 16 875 (± 5955) in 2002 to a peak abundance estimate of 33 488 (± 4455) in 2012. An apparent 20% decline from 2012 to 2021, 33 488 (± 4455) to 26 662 (± 4192), suggests the population abruptly reached carrying capacity due to loss of prey resources. This was particularly evident for humpback whales wintering in Hawai‘i, where, by 2021, estimated abundance had declined by 34% from a peak in 2013, down to abundance levels previously seen in 2006, and contrasted to an absence of decline in Mainland Mexico breeding humpbacks. The strongest marine heatwave recorded globally to date during the 2014–2016 period appeared to have altered the course of species recovery, with enduring effects. Extending this time series will allow humpback whales to serve as an indicator species for the ecosystem in the face of a changing climate.


Introduction
Population monitoring has become increasingly important in conservation biology as anthropogenic activities and climate change continue to impact marine ecosystems [1][2][3][4][5].Abundance trends provide critical insights into the dynamics of animal populations, enabling a better understanding of the ecological interactions and underlying drivers that influence their distribution, abundance and life history [6,7].Abundance estimation is the basis for assessing the current protection status for humpback whales in the USA [8] and Canada [9] and has helped the understanding of anthropogenic effects such as ship strikes on blue whales [10].For marine mammal populations exploited to near extinction, abundance estimation has served as an indispensable tool for decision-makers, supporting the development of effective conservation and management strategies aimed at protecting these animals and the ecosystems they inhabit [11,12].Understanding population trends is essential for addressing the growing challenges faced by marine mammals in a rapidly changing world [13]; however, detecting population changes can be difficult due to the imprecision of most marine mammal abundance estimates and the infrequent intervals at which most populations are surveyed [14].For large cetaceans, such as baleen whales, these issues are largely due to the costs and logistical challenges associated with surveying vast, often remote marine environments, as well as elusive and often cryptic animal behaviour [15].Limited available data may result in poor model parameter estimation and reduced predictive power, ultimately hindering effective conservation and management efforts [16].Furthermore, the infrequent generation of updated population estimates may fail to capture rapid changes in whale populations due to anthropogenic disturbances or environmental fluctuations [17].To address these challenges, we have drawn on an innovative ocean-basin-wide collaboration [18] to cost-effectively generate a time series of relatively royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231462 precise abundance estimates for a long-distance migrating species, the humpback whale (Megaptera novaeangliae) in the North Pacific Ocean.
In the North Pacific, humpback whales are known to breed in waters off Japan, The Philippines, the Mariana Islands, Hawai'i, Mexico and Central America [19][20][21][22].These whales migrate to feed in coastal to continental shelf waters off Russia, the Bering Sea, Alaska, British Columbia, Washington, Oregon and California [20].As with all large cetaceans worldwide, humpback whales in the North Pacific were the target of extensive commercial whaling until late in the twentieth century, with an estimated 31 785 whales taken from 1900 to 1976 [23,24].Humpback whale populations in the North Pacific were severely depleted [24], with the remnant abundance crudely estimated at a low of approximately 1200-1600 individuals [25,26] around the end of humpback whale commercial catches in 1976 [24].The end of commercial whaling created a possibility of full recovery to pre-whaling abundance; it is, however, difficult to define what full recovery means without precise estimates of pre-whaling abundance.Rice [27] suggested that the humpback population was 'stable at about 15 000 [animals]' prior to 1905, and 'was reduced to about 1000' by 1965.However, he gave no supporting evidence for the former figure, and the latter was estimated using catch statistics that were available at the time (either published or provided by the USSR); it is now known that the Soviet statistics were significantly under-reported [28].As Rice stated, no one at the time had attempted to estimate the pristine population size prior to modern whaling, and the number of animals killed in various aboriginal operations prior to the twentieth century is largely unknown.Consequently, there is currently no reliable estimate of prewhaling abundance, although calculating one is a primary aim of an ongoing Comprehensive Assessment by the Scientific Committee of the International Whaling Commission [29].
To understand population status, two previous abundance estimation studies for North Pacific humpback whales were undertaken using mark-recapture (also known as capture-recapture) techniques with synoptic datasets [30,31].These studies depended on visual matching of photographs by skilled technicians to determine marks and recaptures of individual humpback whales based on the unique shape and pigmentation of the ventral surface of the tail (fluke) [32,33].The first photographic markrecapture abundance estimate for humpback whales in the full North Pacific Ocean was based on a historical database of 3650 fluke images gathered from 1990 to 1993 [34].Important western Pacific Ocean feeding areas of Russia, the Aleutian Islands and the Bering Sea were not surveyed during the study, and because it was retrospective, sampling did not follow a systematic design.Using a geographically stratified Darroch method, that study estimated an abundance of approximately 6000 (s.e.= 474) humpback whales in 1992 from three major wintering areas (Mexico, Hawai'i and Japan), but the authors' appraisal of likely biases suggested an adjusted estimate of approximately 8000 [34].Assuming the same coefficient of variation (CV) as the uncorrected estimate, a minimum estimate of the standard error for the bias-corrected estimate is 632.Petersen estimates comparing wintering and summer areas suggest that this may have been an underestimate, but a more accurate Petersen model approach could not be applied uniformly due to the missing data from the unsampled summer areas.Some of the geographical biases in the 1990-1993 sampling were corrected in the 2004-2006 project Structure of Populations, Levels of Abundance, and Status of Humpback whales (SPLASH) [20], when coordinated efforts by over 40 different research groups across the North Pacific obtained fluke identification photographs from all known wintering and summer areas in the North Pacific during three winters and two summers.The experimental design sought to apportion sampling effort in proportion to whale density to representatively sample whales in all areas.The SPLASH project's North Pacific abundance estimate was based on a Chapman-Petersen approach that compared all wintering fluke photographs with all summer fluke photographs [30].In that study, simulations were used to correct biases associated with births and deaths (a common violation of the closed population assumption), not sampling calves and missed matches.Many of these biases were found to be offsetting, and the net estimate of resulting bias was only +3.5%.The bias-corrected SPLASH abundance estimate was 21 063 (CV = 0.04) humpback whales in 2005.
Based on evidence of ongoing recovery indicated by the SPLASH study and other studies, the species was partially delisted from the United States Endangered Species Act (ESA) in 2016 [8], and was reassessed from Threatened to Special Concern by the Committee on the Status of Endangered Wildlife in Canada in 2011 [35].However, multiple current anthropogenic stressors including ship strikes [36], entanglements [37] and climate change [5,[38][39][40] may limit or reverse full recovery.Our 2001-2022 study period included seasons before, during and after periods of a severe marine heatwave [41], with demonstrated effects on regional humpback whale populations [39], and signs of compounding climate change and habitat use conflict such as fisheries interactions leading to exacerbation of these existing threats [42,43].
royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231462 Our current study builds on past photo-identification (photo-ID) efforts mentioned above, incorporating the historical data gathered in those studies, combined with long-term region-specific photo-ID efforts, and making broader use of opportunistic sources, especially whale watch voyage naturalists.The sheer scale of data now available is due to advances in automated image recognition fluke photo-ID matching [44], data management and research collaboration [18].These advances have shifted the abundance estimation paradigm from data scarcity and periodic study to continuous and accessible tracking of the ocean-basinwide population through time.Given that no single research group's data archive or effort could scale and maintain what would be necessary for such an ongoing study, we combined research-based and community science data in a broad collaboration facilitated by the data management Web platform Happywhale.com [18].Where data deficiency previously imposed limits on understanding, this study used a near-real-time population modelling infrastructure and with this, we aimed to evaluate population change with a robust time series of 20 annual abundance estimates from 2002 to 2021.We used the same Chapman-Petersen mark-recapture approach for abundance estimation as was used for the 2004-2006 SPLASH study [30], based on a comparison of three wintering samples with two summer samples for each calculated annual estimation.To address geographically biased sampling in pre-and post-SPLASH years, we used a sub-sampling scheme of the SPLASH years to obtain a bias correction.We examined whether the results of this study might be biased by long-term changes in migratory patterns and in the relatively large proportions of the entire population that winter in Mainland Mexico and Hawai'i, respectively.Finally, we discussed findings of a shift in population trajectories from recovery from industrial whaling to responding to current ecosystem conditions.

Sample collection and matching methods
For this study, all available North Pacific humpback whale photo-ID data from a research collaboration of 46 organizations and 4292 community science contributors were aggregated within a study period of 2001 through 2022, and reconciled in a single dataset through the research collaboration and community science Web platform Happywhale.com.This dataset is fully described by Cheeseman et al. [18].Photo-ID data for the purpose of this study consisted of encounters of individually identified whales, with associated dates and geographical locations.Multiple encounters of the same individual within the same season (summer or winter of a given year) were condensed to one capture occasion (capture), using the observation date closest to the midpoint of the respective season.

Geographical stratification
Each capture was assigned to a geographical stratum based on location, to segregate samples between summer (feeding) areas and winter (breeding) areas and as a method to allocate samples for a bias correction estimation described below.Geographical strata include seven summer areas and six wintering areas (tables 1 and 2; figure 1).Stratification generally followed that used in the SPLASH study [20], with strata based on similarities in migratory destinations.For example, the California and Oregon stratum was separated from the Washington and southern British Columbia stratum because whales from the former migrate almost exclusively to the coasts of Mexico and Central America, whereas a portion of the latter also migrate to the Mexican offshore islands and Hawaiian Island regions.Our large sample size allowed us to discern some migration patterns that were not apparent in the SPLASH study, which resulted in the following changes in the stratification scheme: humpback whales in southern coastal Mexico had similar distribution and migratory destinations to whales in Central America [45,46], so we separated southern Mexico from the rest of Mainland Mexico and pooled it with Central America.The 'Russia' stratum in the SPLASH study included whales from the Kamchatka Peninsula, the Commander Islands and the Gulf of Anadyr.Based on winter migratory destinations, we revised the stratification of Russia as follows: we included the Kamchatka Peninsula as a separate stratum (Kamchatka whales migrated almost entirely to islands in the western Pacific), we combined the Commander Islands with the other western Aleutian Islands to be a new western Bering Sea stratum, we combined northern Russia coastal whales with a stratum that includes the northern and eastern Bering Sea, and samples from the Mariana Islands (which was not sampled in the SPLASH study) were pooled with the western Pacific stratum.Whale distributions in the North Pacific are continuous, with blurred boundaries across most regions (e.g.[47]), with regional delineation used here only to illustrate varied royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231462 sampling over time and to develop bias-correction factors that address non-random sampling, as described below.We do not attempt to estimate abundance for each stratum and our primary abundance estimation method depends on all summer (feeding) areas pooled, and all winter (presumed breeding) areas pooled.Consequently, these geographical stratification changes have minimal impact on actual abundance estimates, serving only to align sample distributions between SPLASH study years and non-SPLASH study years.See electronic supplementary material, I, for stratum boundaries.

Absolute abundance estimation methods: full North Pacific Ocean
Prior to estimating abundance, whale identification data from Happywhale were pre-processed with a custom R script [48].For each capture, sample season was assigned as either summer or winter.As with the SPLASH study, midpoints for each capture period were defined as 1 March (for winter) and 1 August (for summer), hereafter 'season-year'.In the rare cases of an individual encountered in multiple geographical strata in a given season-year, the individual was assigned to only one geographical stratum, based on the closest capture occasion to the seasonal midpoint.
To be comparable to the 2004-2006 SPLASH abundance estimates [30], we used the Chapman biascorrected version of Petersen estimator [49] with wintering area samples ( pooled over three consecutive years) as one capture occasion and with summer area samples ( pooled over the two summers between the three years of winter samples) as the second capture occasion.Modelling has shown that this estimator is robust to sample bias created by heterogeneity in capture probabilities if the factors affecting differential availability are different in the two capture occasions [30].We relied on the different geographical distribution and behaviour of whales between summer and winter to reduce sample bias between mark and recapture samples.Recaptures included those individuals Table 1.Sample size of unique photo-identifications of humpback whale identifications in each of seven summer regions for each season-year.Regions, west to east, include the Kamchatka Peninsula (Kamchatka), the western Bering Sea (WBerSea), the eastern and northern Bering Sea and southern Chukchi Sea (E&NBerSea), the Gulf of Alaska (GulfOfAK), southeast Alaska and northern British Columbia (SEAK&NorthBC), southern British Columbia and Washington (SouthBC&WA), and California and Oregon (CA&OR).Although the simple abundance estimation approach outlined above should work well if the geographical areas were as representatively sampled throughout the time series as during the 2004-2006 SPLASH study, the achieved geographical distribution of samples was uneven in many years (table 1).In particular, the Revillagigedo Archipelago (MxIsl), Baja California, Mexico (MxBaja), Kamchatka, Russia (Kamch) and the eastern and western Bering Sea (BerE and BerW) were intermittently and generally under-sampled.This was expected to bias the abundance estimates outside of the SPLASH years.We therefore used a novel geographical bias correction approach, constructed as follows.

season-year
We assumed that the bias-corrected SPLASH abundance estimate of 21 063 whales [30] was an unbiased estimate for study years (2004)(2005)(2006).We denoted this reference estimate as N SPLASH .Let the photo-ID data underpinning N SPLASH be denoted Y SPLASH .Let N t denote the time series (2002-2021) of abundance values estimated in our analysis and let Y t be the associated photo-ID data.Note, each Y t value represented 3 years of data centred on year t (e.g.Y 2010 are the data from 2009, 2010, 2011).For each t, we subsampled Y SPLASH (call this subsample Y SPLASH, sub(t) ) so that it had the same geographical distribution as Y t (see next paragraph).Y SPLASH, sub(t) was then used to obtain an alternative SPLASH abundance estimate (N SPLASH, alt(t) ).We then calculated the ratios between N SPLASH and each N SPLASH, alt(t) .These ratios (F t = N SPLASH /N SPLASH, alt(t) ) were interpreted as the amount of bias associated with the geographically biased subsamples and were applied as correction factors to the uncorrected annual Chapman estimates (N 0 t ), to obtain corrected estimates, i.e.N t ¼ N 0 t Â F t .To generate subsamples from the SPLASH data (Y SPLASH, sub(t) ) with the same geographical distributions as the Y t , we used a largest common sample size (LCSS) approach.We calculated the .The Chapman-Petersen estimator was used with these data subsets to obtain N SPLASH, alt(t) and N 0 t , which were used to find F t and N t (see previous paragraph).The primary assumptions of this bias-correction approach are that (i) the SPLASH abundance estimate itself is not biased by non-representative geographical sampling, (ii) migration patterns have not changed since the SPLASH period, and (iii) the geographical distribution of whales has not changed during the course of the time series, i.e. that they have remained proportionally distributed throughout the study area as they were during the SPLASH period.The assumption of geographical sampling bias in the SPLASH study was addressed as part of that study [20,30].The second assumption, of migration pattern consistency, was evaluated by comparing migratory destinations of southeast Alaska and northern British Columbia from two samples from different time periods, 2004-2006 and 2019-2021, three seasons each.The third assumption, of consistency of geographical distribution of whales, was evaluated through an ANOVA test of Hawai'i versus Mainland Mexico annual mean relative abundance estimates over the study period to detect shifts in differential population growth rates between the two most sampled wintering areas.

Relative abundance estimation methods: Hawai'i and Mexico
To add resolution to understanding humpback whale abundance over the study period, we estimated a time series of abundances for two wintering areas which were relatively well sampled in most years (Hawai'i and Mainland Mexico, table 2) with consistent research effort over the study period, and which comprised a high proportion of total humpback whale abundance in the North Pacific, with migratory contributions from multiple feeding areas.Chapman-Petersen abundance estimates were made for both wintering areas based on photo-ID samples in consistently well-sampled feeding areas.Feeding areas were selected based on (i) being a substantial migratory contribution to the given wintering area [20] and (ii) having a consistently high sample size of identified whales.The Hawai'i estimate was based on using that wintering area as one sampling occasion and the southeast Alaska/ royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231462 northern British Columbia feeding area stratum as the second sampling occasion.The Mainland Mexico estimate was based on using that wintering area as one sampling occasion and the California/Oregon feeding area stratum as the second sampling occasion.Humpback whales in both of these paired wintering and feeding areas exhibited elements of migratory herd units [50] with wintering areas shared with whales from multiple feeding areas, and feeding areas populated with whales almost entirely from their respective paired wintering areas.As with the ocean-basin-wide model, estimates were based on running 3-year periods with two feeding area samples and three wintering area samples.To interpret these abundance estimates as absolute abundance would have required that the wintering area sample was a random sample of the migratory whales from all feeding areas that contributed to it.Because that assumption was not evaluated in the scope of this study, we treated our estimates as an index of relative abundance.

Standard error estimation
To estimate standard errors (s.e.) in abundance estimates, we used a jackknife approach [51].For both the bias-corrected estimates for the entire North Pacific and for the estimates of relative abundance for the Hawaiian Island and Mainland Mexico strata, jackknife samples were taken by randomly excluding 5% of the fluke identifications (without replacement, after season-year duplicates were excluded for each geographical stratum).The estimation of abundance (including the bias correction) was repeated for each of the 20 jackknife samples, and s.e. were estimated using the standard formula: s:e:(x) where n was the number of jackknife estimates, b x i was the ith jackknife estimate, and x jack was the mean of the jackknife samples.The 2004-2006 SPLASH abundance estimate was used as a reference for correcting the bias for all other sample periods.To include the uncertainty in the SPLASH bias correction estimate (CV = 0.04 [30]), the CV of the bias-corrected abundance was estimated as the square root of the sum of the squared CV of the uncorrected abundance estimate plus the squared CV of the SPLASH estimate.

Sampling distribution
The full 2001-2022 dataset consisted of 30 484 individual humpback whales in 132 684 unique season-year encounters, distilled from 192 869 unique encounters.The number of unique individuals that were photographically sampled each season-year varied from 1940 to 5668, and from 1252 to 3236 for the 21 winters and summers, respectively (table 3).Photo-ID samples were post-stratified into seven feeding areas (table 1) and six wintering areas (table 2).Sample sizes varied among seasons and regions due to varying effort; 20 of 294 season-regions had no effort and zero identifications.Kamchatka and the Mexican Islands regions were consistently under-sampled, with an average of 14 and 49 seasonal identifications, respectively.Hawai'i (833), California/Oregon (648) and southeast Alaska/northern British Columbia (772) averaged the highest seasonal identifications, resulting from both high abundance and high effort.Across all seasons and regions, we found a mean of 448 unique identifications per season, and a maximum of 1798 unique identifications in one season.Sample sizes increased on average fourfold over the years of the study period, with the highest number of unique identifications from 2017 onward (table 3).Further description of this sample distribution is found in [18].

Recapture rates and bias correction
Between-season matches ranged from 4.95% to 11.43% of individuals per season (table 3).Over the full study period, 63% of all individuals were captured in multiple seasons and/or regions, suggesting that a high proportion of the North Pacific humpback whale population was included in the sample [18].While sample sizes were very limited for some regions, no region stood out with a disproportionately low recapture rate.Recapture rates ranged from 56.8% of individuals captured in multiple seasons and/or regions in the Mariana Islands to 87.2% in southern Mexico and Central America.Following the photo-ID methodology of [44], we assumed that 97-99% of potential matches were found.
royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231462 Table 3. Sample size of unique summer and winter area whales (three and two seasons pooled, respectively), number of between-season matches, and resulting North Pacific abundance estimates without and with bias corrections for non-representative geographical sampling.Standard errors are estimated from jackknife samples (n = 20).Between season match rates are (matches between season)/(total unique whales in winter and summer).
ΔYoY is the year-over-year rate of change in abundance, presented as

Abundance estimation and trends
Bias-corrected estimates of humpback whales in the full North Pacific increased from 16 875 (s.e.= 5955) in sample year 2002 to a peak of 33 488 (s.e.= 4455) in 2012 before an inflection with a decline to 26 662 individuals (s.e.= 4192) by sample year 2021 (figure 2 and table 3).Growth from 2002 to 2013 was consistently positive apart from one anomalous year, 2010, averaging 5.9% yr −1 .From 2014 to 2021, estimated abundance consistently declined, averaging −3.0% yr −1 .Hawai'i region relative abundance estimates showed a peak in 2013, then declined to abundance levels equal to the mid-2000s from 2017 onward (table 4 and figure 3).The Hawai'i region growth phase was steeper than the North Pacificwide aggregated dataset, averaging 6.9% yr −1 from 2002 to 2013, then declined at an average of −4.3% yr −1 from 2014 to 2021.By contrast, relative abundance for the Mainland Mexico region grew at an average of 7.1% yr −1 from 2002 to 2015 before appearing to stabilize with an average of 0.9% yr −1 growth from 2016 to 2021 (table 4 and figure 3).

Migration patterns and relative abundance
Analysis of migratory destinations of 2022 individuals found in southeast Alaska and northern British Columbia in 2004-2006 and 2025 individuals in the same area in 2019-2021 found 591 and 698 individuals in wintering areas, respectively, with no significant difference in migratory destinations between the two time periods (χ 2 = 0.263, p = 0.877; table 5).The proportionate relative abundance of whales between Hawai'i and Mexico in 2002-2006 was 78% versus 22%, respectively (table 4 and figure 4).The Hawai'i percentage declined by an average of −0.8% annually over the study period, so that the relative abundances were 65% versus 35% respectively by the 2017-2021 sampling period (ANOVA p < 0.001).

Discussion
This study estimated a total abundance of 26 662 (s.e.= 4192) humpback whales in the North Pacific as of 2021, growing from 2002 to 2021 at an average rate of 3% yr −1 , with a 20% decline from a peak abundance of 33 488   royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231462 in our Hawai'i relative population estimation (figure 3), is consistent with observed changes in oceanographic conditions in feeding areas, specifically the severe northeast Pacific Marine Heatwave (PMH) of 2014-2016 [53], which depressed biological productivity and thus prey availability in northeast Pacific Ocean whale feeding areas.Two other known sources of humpback whale mortality, ship strikes and entanglements, may have contributed to the apparent decline; both have exceeded estimated potential biological removal rates [54,55] and entanglement rates surged on the US West Coast in correlation with the PMH apparently due to habitat compression [43,56,57] but cannot account for an estimated population decline approaching 7000 individuals across the North Pacific in just 9 years (2012-2021).
The PMH, the most extreme in recorded history [41], had profound environmental impacts including elevated sea surface temperatures (by as much as 3 s.d. in some areas with maximum sea surface temperature anomalies sometimes exceeding 3-6°C to depths ranging from 50 to 200 m), a decline in sea surface winds, reduced upwelling and decreases in nutrient-rich water.This led to reduced phytoplankton biomass and restructured zooplankton communities in favour of lower-calorie species [58,59].These changes, in concert with increased metabolically driven food demands of ectothermic forage fish resulting in quality and quantity declines in forage fish [60], and an increased demand for forage fish by large ectothermic predatory fish (e.g.groundfish and salmon), ultimately increased competition beyond the demands already imposed by top predators such as marine mammals and sea birds.The result created what Piatt et al. [61] proposed as an 'ectothermic vise' on the diminished forage resource.High latitude waters within the PMH witnessed a decline in local abundance of humpback whales, where they traditionally feed on krill and small schooling fish such as herring, capelin, sand lance and juvenile salmon [39].Reduced prey biomass of zooplankton combined with increased competition for diminished forage fish resulted in severely diminished presence of calves, prevalence of 'skinny' whales, increased strandings [62], and the absence of individual whales with high site fidelity to specific locales.The impacts were also witnessed in humpback whale breeding grounds including significantly fewer whales and calves of the year (e.g.[38,[63][64][65]).Co-occurring with the negative effects on humpback whales in their feeding grounds, several mass mortalities of fish-eating seabirds were documented suggesting far reaching impacts of the PMH.For example, over 60 000 dead or dying emaciated common murres (Uria aalge) were documented in 2015-2016 from California to Alaska, apparently from starvation, with high rates of reproductive failure [61].The eastern Bearing Sea witnessed shifts in zooplankton community composition and forage fish distribution with a mass mortality of emaciated tufted puffins (Fratercula cirrhata) during the onset of molt, apparently associated with starvation [66].Other species, whose food resources partially overlap with those diminished in association with the PMH, and which also experienced declines in body condition and unusual mortality rates included: various Alaskan large predatory groundfish, California sea lions (Zaolphus californianus) and Guadalupe fur seals (Arctocephalus townsendi) (summarized in [61]).Given all of these co-occurring impacts with the PMH, and because humpback whales feed on both krill and small schooling fish, and may be flexible in their ability to shift their dominant food resource between the two when faced with variations in oceanographic and ecological conditions affecting their food resources [67], they may serve as a valuable bellwether indicator species of ocean-basin ecosystem health.Their lack of recovery by 2021 is notable.
While conservation management has focused on populations defined by breeding areas [8], feeding area resource limitations are likely the greatest determinants of carrying capacity.Regional evidence in Alaska shows variable recovery in feeding areas where local abundance was severely reduced by the PMH, still depressed in Prince William Sound [40] compared to partial recovery in Glacier Bay National Park [39].The documented downward inflection in the population growth curve of North Pacific humpback whales may indicate that population impacts of ecological stressors like the PMH may be broad, and supports evidence of persistence for years after oceanographic conditions return closer to the long-term mean [53].Similar findings of the lingering effects of the PMH have been documented in seabirds in Cook Inlet, Alaska [68].

Bias correction validity and assumptions
The bias-correction approach allowed us to account for variable sampling effort over time and space, necessitating three assumptions.First, the assumption that the SPLASH abundance estimate is unbiased by non-representative geographical sampling is supported by an experimental design for the SPLASH study that pursued representative sampling of all wintering and summer areas known at the time [20,30].The SPLASH study had no effort in the Mariana Islands, which has since been found to be a breeding area for humpback whales [21]; however, after seven seasons of effort there with 37 individuals identified [18], we are confident that this omission from SPLASH sampling effort would royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231462 not significantly change the abundance estimation at the scale of the full North Pacific Ocean.The northwest Hawaiian Islands are now known to be a wintering area [69] with as yet unknown abundance and relationship to the main Hawaiian Islands.Of 37 individuals identified in two northwest Hawaiian Islands voyages of limited scope and duration (2007 and 2019), 24 individuals were resighted in Hawai'i, interchange that suggests a substantial degree of mixing.We believe this is the only North Pacific humpback whale wintering area that may host an under-sampled, site-faithful population, though with a majority of identified individuals also encountered in the main Hawaiian Islands.Further study may thus reveal minimal impact on total abundance estimation statistics.Russian feeding areas were represented by only 102 individuals during SPLASH years [20]; in contrast, this study benefited from data from 2296 individuals [18].This heterogeneity in sampling effort contributed to larger s.e.ranges derived from our jackknife approach in the full North Pacific abundance estimate, compared to narrower ranges in Hawai'i and Mexico relative abundance estimates.
A further consideration of potential sampling bias is of equal detection probability.High accuracy of the AI image recognition tool used for matching [44] meant any fluking whale could be confidently identified irrespective of distinctiveness.Females with calves, however, are known to fluke less often in wintering areas [70] though Craig and Herman demonstrated that extended focal follows of nonfluking individuals would greatly reduce sampling bias [71].The same study found lower wintering resight rates of females.Further, not all humpbacks migrate every year [72].We believe sampling bias is effectively minimized by using multi-year pooled samples for each annual abundance estimation and the methodology of drawing from a dataset of nearly 200 000 identified encounters (192 869 unique encounters distilled to 132 684 unique season-year encounters) that may mimic the high detection rates of extended focal follows [71].
The second assumption, that migratory patterns have not changed, was conclusively found to be valid based on no difference between wintering destinations of over 2000 whales each in time periods near the beginning and end of our study period.This assumption is further supported by demonstrated strong maternal fidelity and natal philopatry in humpback whales with sufficient stability over time to show genetic differentiation [73][74][75].
Validity of the third assumption would require that population trends were uniform between and across all summer and wintering areas during the study period.We found that relative abundance between Hawai'i and Mexico shifted from 78% versus 22% to 65% versus 35% over the course of the study (figure 4) apparently because of different population growth rates (figure 3).This 13% shift in relative abundance between Hawai'i and Mexico, and possible undetected shifts in relative abundance between other study area regions, does violate this assumption and thus may have introduced error in our geographical bias correction factor affecting absolute abundance estimates.We realize that the 13% relative abundance shift is a violation of the assumption of consistent geographical distribution underlying our sampling bias correction methodology; however, we believe its effects on observed abundance trends (as documented by the study) are likely to be minimal because (i) the large samples created by multi-year pooled effort where approximately one-third to half of the total estimated North Pacific humpback population was captured in most 3-year samples and (ii) the bias corrections themselves are relatively small.The overall trends without the bias correction are similar to the biascorrected trends (table 2), although the values are lower in all years except one.We cannot fully discount this as a potential source of bias, particularly after 2013 when the proportion of the population in Hawai'i decreased substantially (figure 4).To address these sources of uncertainty, future study for humpback whales in the North Pacific should include absolute abundance estimation by region.

Relative abundance estimation
We sought to refine the ocean-basin-wide abundance estimate with regional estimates based on the two bestsampled subset wintering areas, Hawai'i and the Mainland Mexico.These subset models could not be considered absolute abundance estimations because each 3-year sample was not assumed to be a closed population or controlled for sample bias such as, in the case of Mainland Mexico, an overlap of Central American whales.These whales are known to co-occur using the region as a migratory corridor [46,76], and therefore the subset abundance estimate sampled a population with undefined boundaries.In Hawai'i, sample effort was concentrated in the Maui-Nui channel and to lesser extents off the islands of Hawai'i, Oahu and Kauai, without a clear definition of whether this is a comprehensive sample of an equally mixed Hawaiian island population.The value of the relative abundance estimations was in describing different population growth trajectories with far lower s.e.estimations than the North Pacificwide model.An apparent sharp decline in humpback whale relative abundance off west Maui (which royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231462 traditionally has hosted the greatest concentrations of humpback whales including calves of the year, e.g.[77][78][79]) was reported after the 2014-2016 PMH, using boat-based transect surveys [38] and passive acoustic monitoring [65], as well as off the Kohala coast of Hawai'i Island using shore-based systematic scans [64].The relative abundance trend measured agrees with data from density estimates by fixed-wing aerial surveys of the entire Maui Nui region for 1993-2003 and 2019-2020 [80], though these surveys did not occur during the period of apparent steep increase then decrease, over the 2005-2016 period.Our results support the conclusions of Frankel et al. [64] that environmental variables likely caused Hawaiian humpback populations to decline due to the PMH.Fluctuating but resilient humpback whale crude birth rates in Hawai'i [38,80] suggest the measured −4.3% yr −1 decline from 2014 to 2021 may be temporary, stabilizing if ocean conditions in feeding areas return to the long-term mean.These trends point to an apparent ecological shift from a half century of recovery response following near-extirpation by industrial whaling to carrying capacity-limited abundance.It is notable that the end of the shift did not appear to occur from population growth but rather from a decrease in carrying capacity triggered by rapid climate change.

Conclusion
The end of the industrial whaling era left oceans largely empty of great whales [23,24].Yet with no species hunted to global extinction during this period, full recovery became a possibility.For humpback whales, steady population growth was the 40-year trend in the North Pacific, and it is to be celebrated that humpback whales appeared sufficiently recovered to qualify for partial delisting from legal protected status regimes in the USA and Canada.Ironically, the timing of legal conservation status changes for the species was coincident with the dramatic ocean-warming-induced decline documented in this study.Some populations of humpback whales may no longer be a priority for endangered species conservation funding [81], but now offer high value as an indicator species of ocean-basin-wide ecosystem health [82] in a world where we can expect increased frequency and severity of marine heatwaves.Past population monitoring efforts offered periodic abundance estimates interspersed with many years with very coarse assessment of population trajectory (figure 2).This study creates a shift to continuous population monitoring with rapid feedback, enabled by an advance in AI-assisted photo-ID [44] and an ocean basin-wide collaborative effort [18].If this costeffective collaborative effort is maintained with funding to gather robust sample sizes from most wintering and summer areas of the species, we can expect to be able to detect shifts in abundance in near real time through ongoing updates of the model with data from each passing year.

Figure 1 .
Figure 1.Locations of photo-identification samples from humpback whales in the North Pacific 2001-2021 colour coded by geographical strata.Summer feeding area strata north of 32°N, west to east are: Kamchatka Peninsula, Russia (black), west Bering Sea (cyan), north and western Bering Sea (green), Gulf of Alaska (yellow), southeast Alaska and northern British Columbia, Canada (blue and purple), southern British Columbia, Canada and Washington, United States (salmon), and California and Oregon, United States (orange).Wintering area strata south of 32°N, west to east are West Pacific ( purple), Hawai'i, United States (blue), Baja California, Mexico (black), Revillagigedo Islands, Mexico (cyan), Mainland Mexico (green), southern Mexico and Central America (yellow).

Figure 2 .Table 4 . 3 .Figure 3 .
Figure 2. Bias-corrected estimates of humpback whale abundance for the full North Pacific Ocean.The red line represents a threeyear moving average of abundance estimates.Vertical lines represent confidence intervals based on ± twice the standard errors (s.e.).S.e. include both the uncertainty in the abundance for the SPLASH years (CV = 0.04, mid-sample year = 2005 [30]) and the uncertainty in the bias-correction process.Inset: North Pacific humpback whale abundance in the context of post-whaling population estimates [25,26] and 1993 estimate [31].royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231462

Figure 4 .
Figure 4. Proportion of Hawai'i and Mainland Mexico relative abundance of humpback whales during study period of 3-year pooled samples 2002 through 2021, with a slope of −0.8% per year.

Table 2 .
Sample size of unique photo-identifications of humpback whale identifications in each of six winter regions for each season-year.Regions, west to east, include the western Pacific (WPac), the Hawaiian Islands (Hawai'i), Baja California, Mexico (MexBajaCal), the Revillagigedo Archipelago, Mexico (MexIsl), Mainland Mexico (MexMld) and Central America and southern Mexico (CenAm&SMex).
number of unique whale identifications in all summer and wintering geographical strata for both the SPLASH period(2004)(2005)(2006)and for the estimation period (t − 1 to t + 1).For each summer and wintering stratum we then calculated the largest sample size that could be taken from both Y SPLASH and Y t .We then subsampled both datasets by randomly sampling this LCSS from each (without replacement) to obtain Y SPLASH, sub(t) and Y t(sub)

Table 5 .
Migratory destinations of 2004-2006 and 2019-2021 southeast Alaska and northern British Columbia humpback whales show no change over the study period.